% ad.m

function [Fz,Fzp,Fnp,Gz,Gzp,Gnp, ...
    Gnpnp,Gznp,Gzpnp,Fzz,Gzz,Gzpz,Gzpzp,Fzpz,Fnpnp,Fznp,Fzpnp,Fzpzp, ...
    Fzpzpzp,Fzpzpnp,Fzpzpz,Fzpzz,Fzpznp,Fzpnpnp,Fzzz,Fzznp,Fznpnp,Fnpnpnp, ...
    Gzpzpzp,Gzpzpnp,Gzpzpz,Gzpzz,Gzpznp,Gzpnpnp,Gzzz,Gzznp,Gznpnp,Gnpnpnp] ...
    =ad(f,g,z_vec,zp_vec,np_vec)

tic

disp(' ');
disp('Constructing analytical derivatives...');

if size(f,1)+size(g,1)~=size(z_vec,2)
    error('Number of equations must equal number of variables.');
end
if length(z_vec)~=size(zp_vec,2)
    error('Variable vectors are constructed inconsistently.');
end

nz = length(z_vec) ;
nn = length(np_vec) ;
nf = length(f) ;
ng = length(g) ;

if length(f)==0 %#ok<ISMT>
    Fz = zeros(0,nz) ;
    Fzp = zeros(0,nz) ;
    Fnp = zeros(0,nn) ;
    Fzz = zeros(0,nz^2) ;
    Fzpzp = zeros(0,nz^2) ;
    Fzpz = zeros(0,nz^2) ;
    Fnpnp = zeros(0,nn^2) ;
    Fznp = zeros(0,nz*nn) ;
    Fzpnp = zeros(0,nz*nn) ;
else
    Fz = jacobian(f,z_vec) ;
    Fzp = jacobian(f,zp_vec) ;
    Fnp = jacobian(f,np_vec) ;
    Fzz = sym(zeros(nf,nz^2)) ;
    Fzpz = sym(zeros(nf,nz^2)) ;
    Fzpzp = sym(zeros(nf,nz^2)) ;
    Fnpnp = sym(zeros(nf,nn^2)) ;
    Fznp = sym(zeros(nf,nn*nz)) ;
    Fzpnp = sym(zeros(nf,nn*nz)) ;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % NOTE ON PARALLEL COMPUTING (PARFOR): everytime when you want to use
    % parfor, you have to open the pools/labs/works (using the command
    % "matlabpool" if you want all the available pools to be opened.
    % Without opening the pool, parfor won't work. Note that you also have
    % to close the pool by using "matlabpool close" before you can open
    % another set of pools. NOTE: for small models, it's best just to use
    % "for" since opening and closing the pools takes time -- for large
    % model, better to use "parfor" (how much better? need to try!!!).
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %matlabpool
    %parfor h=1:nf
    for h = 1:nf
        Fzz(h,:) = transpose(vec(transpose(jacobian(Fz(h,:),z_vec)))) ;
        Fzpz(h,:) = transpose(vec(transpose(jacobian(Fzp(h,:),z_vec)))) ;
        Fzpzp(h,:) = transpose(vec(transpose(jacobian(Fzp(h,:),zp_vec)))) ;
        Fnpnp(h,:) = transpose(vec(transpose(jacobian(Fnp(h,:),np_vec)))) ;
        Fznp(h,:) = transpose(vec(transpose(jacobian(Fz(h,:),np_vec)))) ;
        Fzpnp(h,:) = transpose(vec(transpose(jacobian(Fzp(h,:),np_vec)))) ;
    end
    %matlabpool close
    
    if nargout>18
        % Third order
        Fzpzpzp = sym(zeros(ng,nz^3)) ;
        Fzpzpz = sym(zeros(ng,nz^3)) ;
        Fzpzpnp = sym(zeros(ng,nz^2*nn)) ;
        Fzpzz = sym(zeros(ng,nz^3)) ;
        Fzpznp = sym(zeros(ng,nz^2*nn)) ;
        Fzpnpnp = sym(zeros(ng,nn^2*nz)) ;
        Fzzz = sym(zeros(ng,nz^3)) ;
        Fzznp = sym(zeros(ng,nz^2*nn)) ;
        Fznpnp = sym(zeros(ng,nn^2*nz)) ;
        Fnpnpnp = sym(zeros(ng,nn^3)) ;
        for h = 1:nf
            Fzpzpzp(h,:) = transpose(vec(transpose(jacobian(Fzpzp(h,:),zp_vec)))) ;
            Fzpzpz(h,:) = transpose(vec(transpose(jacobian(Fzpzp(h,:),z_vec)))) ;
            Fzpzpnp(h,:) = transpose(vec(transpose(jacobian(Fzpzp(h,:),np_vec)))) ;
            Fzpzz(h,:) = transpose(vec(transpose(jacobian(Fzpz(h,:),z_vec)))) ;
            Fzpznp(h,:) = transpose(vec(transpose(jacobian(Fzpz(h,:),np_vec)))) ;
            Fzpnpnp(h,:) = transpose(vec(transpose(jacobian(Fzpnp(h,:),np_vec)))) ;
            Fzzz(h,:) = transpose(vec(transpose(jacobian(Fzz(h,:),z_vec)))) ;
            Fzznp(h,:) = transpose(vec(transpose(jacobian(Fzz(h,:),np_vec)))) ;
            Fznpnp(h,:) = transpose(vec(transpose(jacobian(Fznp(h,:),np_vec)))) ;
            Fnpnpnp(h,:) = transpose(vec(transpose(jacobian(Fnpnp(h,:),np_vec)))) ;
        end
    end
end

Gz = jacobian(g,z_vec) ;
Gzp = jacobian(g,zp_vec) ;
Gnp = jacobian(g,np_vec) ;

Gzz = sym(zeros(ng,nz^2)) ;
Gzpzp = sym(zeros(ng,nz^2)) ;
Gzpz = sym(zeros(ng,nz^2)) ;
Gnpnp = sym(zeros(ng,nn^2)) ;
Gznp = sym(zeros(ng,nz*nn)) ;
Gzpnp = sym(zeros(ng,nz*nn)) ;

%matlabpool
%parfor h=1:ng
for h=1:ng
    Gzz(h,:) = transpose(vec(transpose(jacobian(Gz(h,:),z_vec)))) ;
    Gzpz(h,:) = transpose(vec(transpose(jacobian(Gzp(h,:),z_vec)))) ;
    Gzpzp(h,:) = transpose(vec(transpose(jacobian(Gzp(h,:),zp_vec)))) ;
    Gnpnp(h,:) = transpose(vec(transpose(jacobian(Gnp(h,:),np_vec)))) ;
    Gznp(h,:) = transpose(vec(transpose(jacobian(Gz(h,:),np_vec)))) ;
    Gzpnp(h,:) = transpose(vec(transpose(jacobian(Gzp(h,:),np_vec)))) ;
end
%matlabpool close

if nargout>18
    Gzpzpzp = sym(zeros(ng,nz^3)) ;
    Gzpzpz = sym(zeros(ng,nz^3)) ;
    Gzpzpnp = sym(zeros(ng,nz^2*nn)) ;
    Gzpzz = sym(zeros(ng,nz^3)) ;
    Gzpznp = sym(zeros(ng,nz^2*nn)) ;
    Gzpnpnp = sym(zeros(ng,nn^2*nz)) ;
    Gzzz = sym(zeros(ng,nz^3)) ;
    Gzznp = sym(zeros(ng,nz^2*nn)) ;
    Gznpnp = sym(zeros(ng,nn^2*nz)) ;
    Gnpnpnp = sym(zeros(ng,nn^3)) ;
    for h=1:ng
        Gzpzpzp(h,:) = transpose(vec(transpose(jacobian(Gzpzp(h,:),zp_vec)))) ;
        Gzpzpz(h,:) = transpose(vec(transpose(jacobian(Gzpzp(h,:),z_vec)))) ;
        Gzpzpnp(h,:) = transpose(vec(transpose(jacobian(Gzpzp(h,:),np_vec)))) ;
        Gzpzz(h,:) = transpose(vec(transpose(jacobian(Gzpz(h,:),z_vec)))) ;
        Gzpznp(h,:) = transpose(vec(transpose(jacobian(Gzpz(h,:),np_vec)))) ;
        Gzpnpnp(h,:) = transpose(vec(transpose(jacobian(Gzpnp(h,:),np_vec)))) ; 
        Gzzz(h,:) = transpose(vec(transpose(jacobian(Gzz(h,:),z_vec)))) ;
        Gzznp(h,:) = transpose(vec(transpose(jacobian(Gzz(h,:),np_vec)))) ;
        Gznpnp(h,:) = transpose(vec(transpose(jacobian(Gznp(h,:),np_vec)))) ;
        Gnpnpnp(h,:) = transpose(vec(transpose(jacobian(Gnpnp(h,:),np_vec)))) ;
    end
end

disp(toc);

end






